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ABSTRACT 

Gravitational wave signals from coalescing Massive Black Hole (MBH) binaries could be used as 
standard sirens to measure cosmological parameters. The future space based gravitational wave 
observatory Laser Interferometer Space Antenna {LISA) will detect up to a hundred of those events, 
providing very accurate measurements of their luminosity distances. To constrain the cosmological 
parameters we also need to measure the redshift of the galaxy (or cluster of galaxies) hosting the 
merger. This requires the identification of a distinctive electromagnetic event associated to the binary 
coalescence. However, putative electromagnetic signatures may be too weak to be observed. Instead, 
we study here the possibility of constraining the cosmological parameters by enforcing statistical 
consistency between all the possible hosts detected within the measurement error box of a few dozen 
of low redshift {z < 3) events. We construct MBH populations using merger tree realizations of the 
dark matter hierarchy in a ACDM Universe, and we use data from the Millennium simulation to 
model the galaxy distribution in the LISA error box. We show that, assuming that all the other 
cosmological parameters are known, the parameter w describing the dark energy equation of state can 
be constrained to a 4-8% level {2a error), competitive with current uncertainties obtained by type la 
supernovae measurements, providing an independent test of our cosmological model. 
Subject headings: black hole physics - gravitational waves - cosmology: cosmological parameters — 
galaxies: distances and redshifts - methods: statistical 



1. INTRODUCTION 

The Laser Interferometer Space Ant enna ( LISA , 
iDanzmann fc the LISA Study TeanJ |1997[ ) IS a space 
based gravitational wave (GW) observatory which is ex- 
pected to be launched in 2022+ . One of its central scien- 
tific goals is to provide information about the cosmic evo- 
lution of massive black holes (MBHs). It is, infact, now 
widely recognized that MBHs are fundamental building 
blocks in the process of galaxy formation and evolution; 
they are ubiquitous in nearby galaxy nuclei (see, e.g., 
iMagorrian et al.|[l998f ). and their masse s tightly corre- 
late with the properties of their host (jGiiltekin et al.l 
120091 and references therein). In popular ACDM cos- 
mologies, s tructure formation p roceeds in a hierarchi- 
cal fashion () White fc Reed [1978( 1 . through a sequence of 
merging events. If MBHs are common in galaxy cen- 
ters at all epochs, as implied by the notion that galaxies 
harbor active n uclei for a short period of their lifetime 
(jHaehnelt fc Rc es 1993), then a large number of MBH bi- 
naries are expected to form during cosmic history. LISA 
is expected to observe the GW driven inspiral and fi- 
nal coalescence of such MBH binaries out to very high 
redshift with high signal-to- noise ratio (SNR), allowing 
very accurate measurements of the binary parameters. 
The collective properties of the set of the observed coa- 
lescing binaries will carry invaluable information for as- 
trophysics, making pos sible to constrain models of MBH 
formation and growth ([Plowman et al.|[20Tot iGair et all 
[2OT0t ISesana et al.l[20m 

Besides astrophysical applications, coalescing 
MBHs could be used as standard sirens (jSchutd 
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of the GW signals allows us to measure the luminosity 
distance with a precision of less than a percent at 
redshift z = \ (neglecting weak lensing). However, 
we need an electromagnetic identification of the host 
in order to measure the source redshift and be able 
to do cosmography. If the event is nearby {z < 0.4), 
then we have a very good localization of the source on 
the sky and we can identify a single cluster of galaxies 
hosting the merger. As we go to higher redshifts, LISA 
sky localization abilities become quite poor: a typical 
sky resolution for an equal mass 1O^M0 inspiralling 
MBH bina ry at z = 1 is 20-30 a rcminutes a side at 2a 
jTrias & Si'Sei[200l iLang fc Hug hes 2001 lArun et al.l 
l2009al ). which is in general not sufficient to uniquely 
identify the host of the GW event. There is, therefore, 
a growing interest in identifying putative electromag- 
netic signatures associated to the MBH binary before 
and/or afte r the final GW driven coalescence (for a 
review, see iSchnittmanI 120101 . and references therein). 
Electromagnetic anomalies observed before or after 
the coalescence within the LISA measurements error 
box may allow us to identify the host and to make a 
redshift measurement. However, most of the proposed 
electromagnetic counterparts are rather weak (below the 
Eddington limit), and in case of dry mergers (no cold 
gas efficiently tunneled into the remnant nucleus) we 
do not expect any distinctive electromag netic transient . 
This brings us back to the original idea bv lSchut3 (|1986[ ) 
to consider each galaxy within the LISA measurement 
error box as a potential host candidate. The idea is 
that, by cross-correlating several GW events, only one 
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galaxy (cluster of galaxies) in each error box will give us 
a consistent set of parameters describing the Universe. 
Th e effectiveness of thi s method has been demonstrated 
by iMacLeod fc Hoganl ()2008l ) in the context of the 
Hubble constant determination by means of low redshift 
{z < 0.2) extreme mass ratio inspirals. 

We use the hierarchical MBH form ation model sug- 
gested bv lVolonteri &: BegelmanI (|201Cl[ ) to generate cata- 
logs of coalescing MBH binaries along the cosmic history. 
This model predicts ^ 100 MBHs mergers observable by 
LISA in three years, in the redshift range [0 : 5]. We 
do not use sources beyond redshift z = 3 due to difficul- 
ties of measuring galaxy redshifts beyond that thresh- 
olcfl. We model the galaxy distrib ution in the Univers e 
using the Millennium simulation (iSpringel et al.|[2005[ ). 
For each coalescing MBH in our catalog, we select a host 
galaxy in the Millennium run snapshot closest in red- 
shift to the actual redshift of the event. For each galaxy 
in the snapshot, we compute the apparent magnitude in 
some observable band, and we create a catalog of red- 
shift measurements of all the observable potential host 
candidates. Note that typical observed mergers involve 
10"^ - 10^ Mq MBHs, which im plies (using the black hole 
mass-bulge relations, see, e.g., IGiiltekin et"an[200l rel- 
atively light galaxies. However, observed galaxies are 
heavy due to selection effects: roughly speaking, mass 
reflects luminosity, so that at high redshifts we can ob- 
serve only very massive (luminous) galaxies. Therefore, 
the actual host might not be (and often is not) among 
the observed galaxies. The important fact is the self- 
similarity of the density distribution: the local density 
distribution for all galaxies and the density distribution 
for heavy galaxies are quite similar, which allow us to 
infer the likelihood of the host redshift on the basis of 
redshift measurements of the luminous galaxies only. 

We assume that the GW source parameter measure- 
ments (GW likelihoods) are represented by multivari- 
ate Gaussian distributions around the true values, with 
the variance-covariance matrix defined by the inverse of 
the Fisher matrix. This is a good approximation in the 
case of Gaussian instrumental noise and large SNR. At 
z > 0.25 the uncertainty in the luminosity distance (-Dl) 
is dominated by weak lensing due to the extended dis- 
tribution of dark matter halos between us and the GW 
source. In this paper we combine the luminosity distance 
errors given by GW measurements and weak lensing, re- 
ferring to them as GW-fWL errors. We use two esti- 
matio ns of the weak lensing erro r (i) fr om lShapiro et ahl 
((20T0h and (ii) from lWang eFall (l200l) . 

In order to evaluate the error box we need to assume 
some prior on the cosmological parameters. In this ex- 
ploratory study, we assume that we know all the cosmo- 
logical parameters but the effective equation of state for 
the dark energy, described by the parameter w (which 
could be the case by the time LISA will fly). In a fol- 
low up paper we will relax this assumption by including 
also the Hubble constant and the matter and dark energy 
content of the Universe as free parameters. We take the 
prior range for w from the seven-year WMAP analysis 
(jKomatsu et al.|[2010D . We show that using statistical 
methods w can be constrained to a 4-8% level {2a error 

^ There are other reasons for not going beyond z = 3 which we 
will discuss later. 



), providing an effective method for estimating the dark 
energy equation of state. We also show that this result 
depends weakly on the prior range and could serve as an 
independent way of measuring the dark energy equation 
of state, with respect to canonical methods employ ing 
observations of type la supernovae (jRiess et al.lll998[) . 

The paper is structured as follows. In Section [5] we 
spell out explicitly all the details of the adopted cosmo- 
logical model and of the Bayesian analytical framework. 
In Section |3] we give more insights on the MBH popu- 
lation model and on the galaxy distributions extracted 
from the Millennium database. In Section!?] we describe 
our simulated GW and electromagnetic observations. We 
give results of our simulations under different assump- 
tions about weak lensing, depth of the follow up electro- 
magnetic surveys, etc. in Section [5j We summarize our 
findings in Section [51 

2. ANALYTICAL FRAMEWORK 

2.1. Cosmological description of the Universe 

We assume the standard ACDM model, which de- 
scribes our Universe as the sum of two non-interacting 
components: (i) a pressureless component correspond- 
ing to all visible and dark matter, (ii) a dark energy 
component with current effective equation of state cor- 
responding to the A— term p = — e. Current estimates 
based on SNla observations and anisotrop y measure- 
ments in the cosmic micro wave background (jRiess et al.l 
119981: IKomatsu et"aIll2010D tell us that about 70% of the 
Universe energy content is in the form of the dark energy. 
The evolution of the Universe is therefore described by 
the expansion equation 



r!t'„(l + ^)' + f^deCxp 3 



dz 



1+Lj{z) 
l + Z 



(1) 

where H = a/a {a being the lengthscale of the Universe) 
is the Hubble expansion parameter and is its cur- 
rent value (t ~ 0), il,n and fi^je are the ratios of the 
matter density and the dark energy density to the criti- 
cal density, and uj{z) describes the effective dark energy 
equation of state as a function of z. We assume that 
the Universe is spatially flat, the luminosity distance is 
therefore computed as 



(1 + ^) 



dz' 



(2) 



In our simulations we fix all parameters (assuming that 
they are known exactly) to the currently estimated mean 

values: Hq = 73.0 km x s"^ x Mpc~\ Vim = 0.25, fide = 
0.75. We also simplify the form of uj{z) for which we will 
assume uj — —1 ^ w, where is a constant 0. We choose 
the value it; = to simulate our Universe which is what 
has been used in the Millennium simulation (see below). 

2.2. Methodology and working plan 

Our aim is to show that we can constrain w via GW 
observations of spinning MBH binaries, using a Bayesian 
framework. Let us consider j = l,..,iVcv GW obser- 
vations. For each event we can infer the probability of 

^ Here we use notations for the d ark energy equation o f state 
adopted in the WMAP data analysis IIKomatsu et al.ir2010l) . 
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a parameter w, given the collected data s, using Bayes 
theorem: 

P,H.) = M!^:1^. (3) 

Here Pj (w\s) is the posterior probability of the parameter 
w, Pj{s\w) is the likelihood of the observation s given the 
parameter w, po{w) in the prior knowledge of w and Ej 
is defined as 



Po{w)Pj (s\w)dw. 



(4) 



The likelihood Pj{w\s) must be appropriately specialized 
to our problem. We want to exploit GW observations 
to constrain w through the distance - redshift (Dl ^ z) 
relation as given by 

• The distance is provided by the GW observa- 
tions: the GW signal carries information about the 
parameters of the binary, including its location on 
the sky and its luminosity distance. All those pa- 
rameters c an be extracted using l atest data analysis 
meth ods (jPetiteau et al.l 120101 : ICornish fc Porteii 
120071 ). The measurements errors are encoded in 
the GW likelihood function C{Dl, 0, (j), A), where 
{9, ^} are the ecliptic coordinates of the source and 
A represents all the other parameters characterizing 
MBH binary (spins and their orientation, masses, 
orientation of the orbit and MBHs position at the 
beginning of observations). When estimating 
weak lensing can not be neglected. In fact the error 
coming from the weak lensing (causing fluctuations 
in the brightness of the GW source which gives an 
uncertainty in the luminosity distance) dominates 
over the GW error starting from redshift z ~ 0.25 
(see figure 121) • 

• The redshift measurement does not rely on any 
distinctive electromagnetic signature related to the 
GW event. We extract a redshift probability dis- 
tribution of the host from the clustering properties 
of the galaxies falling withing the GW+WL error 
box. This defines an astrophysical prior piO, (j), z) 
for a given galaxy in the measurement error box 
to be the host of coalescing binary. To translate 
the measured D l and uncertainty A_D l of the GW 
event into a corresponding z and Az for the can- 
didate host galaxies in the sky we use the prior 
knowledge of po{w) obtained from WMAP. 

The likelihood in equation ([3]) can therefore be written 

as 

Pjislw) = / Cj loLiz, w), 0, (p, X] p{X)pjiO, (f), z) dX dO dq 



where we have introduced the priors p{X) on the parame- 
ters A (which we assume in this paper to be uniform). It 
is convenient to change the variable of integration from 

Through the paper, with GW likelihood wc mean the likelihood 
of the LISA data to contain the GW signal with a given parameters, 
not to be confused with the likelihood Pj{s\w) defined in the Bayes 
theorem. 



Dl to z. Since we have assumed uniform priors on A, we 
can marginalize the likelihood over those parameters to 
obtain: 



P,{s\w) 



TTj [Di^{z, w),9, (/)] Pj{0, (p, z) dO d(/) dz, (6) 



where we denoted the marginalized GW likelihood as 
TTj [£'i(z, w), 0, (/)]. Practically, we limit the integration 
to the size of the error box (in principle the integration 
should be taken over the whole range of parameters but 
we found that considering the 2a error box is sufficient). 

We assume that the error in luminosity distance from 
the weak lensing is not correlated with the GW mea- 
surements, hence the integral in equation © can be per- 
formed over the sky {{0, 0}) first, and then over the red- 
shift. We also found that the correlation between Dl 
and the sky position coming from the GW observations 
is not important for events at z < 0.5. Plugging equation 
^ into equation (O defines the posterior distribution of 
w for a single GW event (as indicated by the index j). 
Assuming that all A^ov GW events are independent, the 
combined posterior probability is 



P{w) 



(7) 



IPo{w)Ylj=iP3{s\w)dw 
To evaluate w through equation ([7]) we therefore need: 

• a MBH binary population model defining the prop- 
erties of the A'ov coalescing systems; 

• the spatial distribution of galaxies within a volume 
comparable with the combined GW+WL measure- 
ment error box; 

• the measurement errors associated to GW ob- 
servations of coalescing MBH binaries (defining 

C,iDL,e,(f>,X)); 

• an estimation of spectroscopic survey capabilities 
to construct the galaxy redshift distribution within 
the GW-I-WL measurement error box (defining 

We will consider these points individually in the next two 
sections. 

3. ASTROPHYSICAL BACKGROUND 

3.1. Massive black hole binary population 

To generate populations of MBH binaries in the Uni- 
verse, we use the re sults of merger tree s imulations de- 
scribed in details in IVolonteri et all (|2003[ ). MBHs grow 
hierarchically, starting from a distribution of seed black 
(^2, holes at high redshift, through a sequence of merger and 
accretion episodes. Two distinctive type of seeds have 
(5^een proposed in the literature. Light (M ^ IOOMq) 
seed are thought to be the rernnant of Population HI 
(POPIII) stars (jMadau fc ReesI 1200 It) , whereas heavy 
seeds form following instabilities occurring in massive 
protogalactic disks. In the model proposed by Begel- 
man Volonteri & Rees (jBegelman et al.l 120061 hereafter 

^ Here this corresponds to the projection of the Fisher matrix 
to three dimensional parameter space of sky location 8, <j> and lu- 
minosity distance ■ 
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BVR model), a 'quasistar' forms at the eenter of the 
protogalaxy, eventuahy cohapsing into a seed BH that 
efficiently accretes from the quasistar envelope, result- 
ing in a final mass M ^ few xlO*M0. Here we use 
the niodel recently suggested by Volontcri & Begcl- 
man (jVolonteri fc Begelmanl[2"010i . hereafter VB model), 
which combines the two above prescriptions by mix- 
ing light and heavy initial seeds. This model predicts 
~ 30 — 50 events per year in the redshift range < z < 3, 
relevant to this study. The dashed blue lines in figure [1] 
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Fig. 1. — Population of coalescing MBH binaries in three years. 
Top left panel: total redshifted mass distribution; top right panel: 
mass ratio distribution; lower panel: redshift distribution. Color 
and linestyle codes are labeled in the figure. 

show the redshifted total mass (M^, = (Mi -I- M2)(l 4- z), 
being M\ > AI2 the rcstframe masses of the two MBHs. 
upper-left panel), mass ratio (q = M2/M1, upper-right 
panel) and redshift (lower panel) distribution of the MBH 
binaries coalescing in three years, as seen from the Earth. 
The model predicts ^ 40 coalescences in the redshifted 
mass range IO^Mq < < 10"^ Mq, almost uniformly 
distributed in the mass ratio range 0.1 < g < 1, with a 
long tail extending to g < 10"'^. For comparison we also 
show the population expected by a model featuring heavy 
seed only (BVR model, green dotted-dashed lines), and 
by an alternative VB type model (labeled VB-opt for op- 
timistic, red long-dashed li nes) with a boo sted cfficicnc; 
of heavy seed formation (see lVolonteri fc B cgclman 201 
for details). It is worth mentioning that these models 
successfully reproduce several properties of the observed 
Universe, such as the present day mass density of nuclear 
MBHs and the optical and X-ray luminosity function s 
of quasars (|Malbon et al.l 120071: iSalvaterra et al.ll2007[ ). 
The BVR and the VB-opt models predict MBH popula- 
tion observables bracketing the current range of allowed 
values. The VB-opt model, in particular, is borderline 
with current observational constraints on the unresolved 
X-ray background, and it is shown here only for compar- 
ison. In the following, we considered the VB model only, 
which fits all the relevant observables by standing on the 



conservative side. 

We performed 100 Monte Carlo realizations of the 
population of MBH binaries coalescing in three years. 
Each realization takes into account the distribution of 
the number of events and MBH masses with the redshift 
as predicted by the VB model. Other parameters (like 
time of coalescence, spins, initial orbital configuration) 
are chosen randomly using uniform priors over the ap- 
propriate allowed ranges. 

3.2. Galaxy distribution 

To simulate the galaxy distribution in the Universe 
we use the data produced by the Virgo Consortium 
publicly available at http: / /ww w. g- vo.org/Millennium, 
These data are the result of the implementation of semi- 
analytic models for galaxy formation and evolution into 
the dark matter (DM) halo r nerger hierarchy gen erated 
by the Millennium simulation (jSpringel et al.ll20 05f). The 
Millennium run is a N-body simulation of the growth of 
DM structures in the expanding Universe starting from 
a Gaussian spectrum of initial perturbations in the DM 
field at high redshift, which successfully reproduced the 
net-like structure currently observed in the local Uni- 
verse. The simulation has a side-length of « 700 Mpc 
(co-moving distance), and its outcome is stored in 63 
snapshots evenly separated in log(z), enclosing all the 
properties of the DM structure at that particular time. 
Semi-analytical models for galaxy formation are imple- 
mented a posteriori within the DM structures predicted 
by the simulation. Such models have been successful 
in reproducing several observed properties of the local 
and the high redshift Un i verse (see , e.g., iBower et al.l 
[200I IDe Lucia fc BlaizotI [2001 . Here we use the im- 

S ilementation perfor med by Bertone and collaborators 
Bertone et al. | [20071 ) . which is a refinement of the origi- 



nal implementation bv lDe Lucia fc Blaizot (2 0(37[ ). 

For each coalescing MBH binary, we choose the snap- 
shot closest in redshift. Within the snapshot we choose 
the host of the GW signal according to a probability 
proportional to the number density of neighbor galax- 
ies rigai. Such assumption comes from the fact that two 
galaxies are needed in order to form a MBH binary, and 
we consider that the probability that a certain galaxy was 
involved in a galaxy merger is proportional to the num- 
ber of neighbor galaxies. We consider to be neighbors of 
a specific galaxy all the N{R) galaxies falling within a 
distance 



R = aTniz), 



(8) 



where cr = 500 km is the typical velocity dispersion of 
galaxies with respect to the expanding Hubble fiow, and 
Th{z) is the Hubble time at the event redshift. The num- 
ber density of neighbor galaxies is then simply written 
as rigai = 3N{R)/{4:TrR^). When we choose the merger 
host, we compute rigai considering all the neighbor galax- 
ies, without imposing any kind of mass or luminosity se- 
lection. In this case rigai = 7^total• However, when we will 
construct the probability of a given observable galaxy to 
be the host of the merger (i.e. the astrophysical prior 
Pj{0,4>, z)), we will have to compute rigai according to 
the number of observed neighbors, because this is the 
only thing we can do in practice when we deal with an 
observed sample of galaxies (see Section |4?2|) . 
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4. SIMULATING THE OBSERVATIONS 

4.1. Gravitational wave observations: shaping the error 

box 



As we mentioned in Section [3. 1[ we drawn hundred re- 
alizations of the MBH binary population from the VB 
model. Each realization contains 30 to 50 events in the 
redshift range [0 : 3]. The total mass, mass ratio and 
rcdshift distributions of the events are shown in the fig- 
ure [TJ In order to simulate GW observations, the binary 
sky location is randomly chosen according to a uniform 
distribution on the celestial sphere, the coalescence time 
is chosen randomly within the three years of LISA op- 
eration (we assume 3 years as default mission lifetime), 
the spin magnitudes are uniformly chosen in the interval 
[0 : 1] in units of mass square, and the initial orienta- 
tions of the spins and of the orbital angular momentum 
are chosen to be uniform on the sphere. More detailed 
descriptio n of the model for GW signal used in this paper 
IS given m iPetiteau et all (|2010( ). 

The GW likelihood £ needed in equation ([5]) is ap- 
proximated as a multivariate Gaussian distribution with 
inverse correlation matrix given by the Fisher informa- 
tion matrix (FIM) : 



'{s-h\s-h) 



(9) 



Here 9' is the vector of the parameters characterizing the 
spinning MBH binary, 0* are the maximum likelihood es- 
timators for those parameters which are assumed to cor- 
respond to the true values (no bias), and Fy = {h^i\hj) 
is the FIM, where the commas correspond to derivatives 
with respect to the parameters. This is a reasonable 
approximation due to the large SNR (for more details 
on the FIM and its applicability see IVallisneril 120081 ) . 
Our unce rtainties on estimated pa rameters are consis - 
tent with iLang fc Hughel (|2009l ). iBabak et all (|20lol ) 
and IPetiteau et al.l ( 20101) . We did not include higher 
harmonics (only the dominant, twice the orbital fre- 
quency) as they only slightly improve parameter esti- 
mation for precessing binaries. However including higher 
harmonics in the GW signal model is important in case of 
the small spins and low precession (when spin s are almost 

i anti) aligned with the orbital momentum, iLang et al.l 
Mil)). We use truncated waveforms corresponding to 
the inspiral only. However the addition of merger and 
ring-down will furth er reduce the locahzatio n error due 
to the higher SNR (jMcWilliams et al.ll201Cil ). This er- 
ror is usually an ellipse on the sky but we simplify it by 
choosing the circle with the same area. 

For the luminosity distance measurement we need to 
take into account the weak lensing. We assume the 
we ak lensing erro r to b e Gaussian with a a given by 
(i) iShapiro et al.l (|2010( ). Such assumption is rather 
pessimistic; we als o tried the prescription given by (ii) 
iWang et alj (|200^ . which gives smaller errors, but still 
larger than the level that may be a chieved after mitiga - 
tion through shear and flexion maps (jHilbert et al.ll2?nol ). 
Both of those estimations are represented in figure |2| as 
(i) dark (red online) circles and (ii) light (orange online) 
squares correspondingly. The median error in due to 
GW measurements only is given by the solid black line. 
The combined error for model (i) is given by the upper 
(blue) circle-line curve, and for model (ii) by the lower 



(green) square-line curve. We consider our setup to be 

0.06 I ^ — 



GW error 

Weak lensing en^or #1 
Weak lensing error #2 
' Combined error #1 
I Combined error #2 



J 0.03 




Fig. 2. — Relative error in the lumin osity distance due to 
weak lensing from (i) IShapiro et al.l 1 12010 ) (circles) and from (ii) 
I Wane et al.l I I2002I) (squares). The black solid line is the median 
error due to GW measurements only; the solid-circle and the solid- 
square lines are for the combined errors under assumptions (i) and 
(ii) respectively (see text). 

conservative in the estimation of the weak lensing effects. 
The main aim of this work is to build a reasonable setup 
for what could be observed by the time LISA will fly, 
and make a first order estimation of LISA capabilities 
to constrain the dark energy equation of state. We will 
address non-Gaussianity of the weak lensing as well as 
other corrections to the model to make it more realistic 
in a follow up paper. 

We consider an error box size corresponding to 2a of 
the measurement errors in the sky location (csky) and in 
the source distance as evaluated by the FIM plus weak 
lensing uncertainties. For observational purposes, the 
dimensions of this error box are AfJ = 2(Tsky and Az. 
For the latter we also include the uncertainty given in the 
Di — z conversion due to the error (prior) on w, po{w). 

Let us summarize how we construct an error box in 
practice, as, for example, the one illustrated in figure 

• We select the closest Millennium snapshot to the 
event in redshift. 

• We pick a galaxy (red dot) in the snapshot with a 
probability given by the local galaxy number den- 
sity ntotal- 

• We construct around the galaxy an error box given 
by Ari and Azqw+wl, and the galaxy can lie any- 
where with respect to this error box (blue cylinder). 

• We expand the error box along the direction of the 
observer both sides by Az given by the uncertainty 
in w (green cylinder). 

• According to some prescription,which we will de- 
scribe in the next section, we select observable 
galaxies in the error box (brown dots). 

As shown in figure [31 we interpret one of the direc- 
tions in the Millennium snapshot as distance from the 
observer, and convert the comoving distance in rcdshift. 
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Fig. 3. — Example of error box (cylinder) in part of the Millen- 
nium snapshot (cube with unit in Mpc). The blue cylinder is the 
measurement error box and the green one also considers the prior 
on w. The black big dot is the host and the brown small dots are 
the selected galaxy candidates. 

We assume a periodic expansion of the Millennium data 
in order to fit large error boxes. Note that the original 
Millennium simulation also assumes the same periodic- 
ity in the distribution of the matter. The size of the 
error box at high redshift covers a significant fraction 
of the simulation box so we do not go beyond the red- 
shift z = 3 (as we will show later, spectroscopic observa- 
tions at such high redshifts will be impractical anyway). 
Together with larger error boxes, we have a nonlinear 
increase in the number of events at high redshift. To re- 
duce the overlap between error boxes corresponding to 
different GW events we choose cylinders with random 
orientations. 

Figure |4] shows an example of the resulting weighted 
distribution of galaxy redshifts (with weight proportional 
to the local density ritotai)- It is a projection of the 
dumpiness along the line of sight which is also propor- 
tional to the probability distribution of z for the event. 
The probability distribution of w for the event will be 
directly related to this result. We noticed that there is a 
very large number of underdense regions and several very 
dense superclusters. The probability of a galaxy with a 
low local density to host a merger is very low but there 
is a huge number of such galaxies, and we found that the 
probability of the host to be in (super) clusters is similar 
to that of being in a low density region. As we will see 
later in the result section, this may cause a very wrong 
estimation of w for some individual GW event. 

4.2. Redshift measurements through spectroscopic 
surveys 

To get a statistical measurement of w we need to ex- 
ploit the clustering of the galaxies falling within the error 
box (which defines the astrophysical prior pj(9^(f),z) in 
equation ([S])). It is therefore necessary to get efficient 
redshift measurements of thousands of galaxies within a 
small field of view (FOV) : the information we seek is en- 
closed in the redshift distribution of such galaxies. We 
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z 

Fig. 4. — Distribution of the weighted galaxies with the redshift. 
The green dashed vertical line is the redshift of the host galaxy. 



stress here that we are not looking for a distinctive elec- 
tromagnetic counterpart to the GW event. In fact, the ac- 
tual host of the coalescing binary may not even be observ- 
able. Typical masses of our binarie s are 10^ — lO^Mfr-,. Us - 
ing MBH-bulge scaling relations (jGiiltekin et al.l[2009( ). 
such MBHs are expected to be hosted in galaxies with 
stellar mass 10® — lO^^M©, i.e., in dark matter halos 
with total mass < 10^^ Mq. The Millennium run mass 
resolution is ^ IO^Mq, meaning that typical host struc- 
tures are formed by less than 100 particles. Unfortu- 
nately, the Millennium run is severely incomplete in the 
expected mass range of LISA MBH binary hosts. Here 
we do not attempt to exploit any MBH-host relation to 
select the host of our GW event; the probability of be- 
ing a host is only related to the local number density 
of neighbor galaxies ntotai- Such assumption relics on 
the concept of self- similarity of the galaxy clustering at 
different mass scales: typical LISA MBH binary hosts 
cluster in the same way as more massive galaxies. We 
checked this assumption by comparing the spatial distri- 
bution of galaxies in different mass ranges {10^ — IQ^'^Mq, 
10^° - 1O"M0, 10" - within simulation snap- 

shots at different redshift, and we postulate that this self- 
similarity extends to lower masses, below the Millennium 
run resolution. This point is crucial for two reasons: (i) 
especially at z > 1, we will be able to get only spec- 
tra of luminous (massive) galaxies, and we need to be 
confident that their spatial distribution mimics that of 
lighter galaxies that may host the GW event but are ob- 
servable in the spectroscopic survey; (ii) the number of 
observable galaxies in the error box may be too large 
anyway (> 10*) to efficiently complete a spectroscopic 
survey on the full sample: self-similarity allows us to get 
the clustering information we need by getting spectra of 
the brightest objects only. 

At z = 1, the typical number of galaxies enclosed in 
the 2a error box described above is in the range 10''^ — 10^. 
However, not all of them arc bright enough to get use- 
ful spectra. The se mianalytic galaxy evolution model 
(jBertone et al.l[2"007f l implemented on top of the Millen- 
nium run returns the stellar mass of each galaxy, and 
the absolute bolometric magnitude Mb. By knowing 
the redshift, and by using standard galactic templates 
one can therefore compute the apparent magnitude in 
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a given band, by assu ming the appropriate k correction 
(|Oke fc SandagelfT96l . Here we use the R band appar- 
ent magnitu de fo r ilhistrative purposes, and we adopt 
the relation (jZombc ck 1990) 

Mb ^ -51og(zc/i7o) - 1-086Z - 25 + m,. + 0.6, (10) 

where 0.6 is the k correction. For each galaxy we com- 
pute rrir and we simulate spectroscopic surveys at differ- 
ent thresholds = 24, 25, 26. We stress here that the 
GW host was chosen among all the galaxies falling in the 
error box, and therefore may not (and usually does not) 
belong to the observed sample. We then assume that for 
each galaxy satisfying the survey threshold we get an ex- 
act spectroscopic redshift, and we combine the redshift 
distribution of several error boxes to get a statistical es- 
timation of w. In practice, each redshift estimation will 
come with a measurement error, and an intrinsic error 
due to the proper motion of the source with respect to 
the Hubble flow. Both errors are however of the order of 
Az/z < 10~^, well below the typical redshift scale corre- 
sponding to spatial clustering of structures (Az ~ 0.01, 
see figure S]) we need to resolve. 

Our method does not rely on the observation of a 
prompt transient associated to the MBH binary coa- 
lescence to identify the host galaxy. Nevertheless, get- 
ting thousands (or tens of thousands) of spectra in a 
small field of view requires a dedicated observational pro- 
gram. Thanks to multi-slit spectrog raphs such as VI- 
MOS at VLT (ILe Fevre et al.l [2001 and DEIMOS at 
Keck (iFaber et al.ll2003l ). fast deep spectroscopic surveys 
of relatively large FOV are now possible . For example, 
the o ngoing VIMOS VLT deep survey (iLe Fevre et al.l 
|2005[) . took spectra of > 10000 galaxies, mostly in the 
redshift range < z < 1.5, within a FOV of 0.61deg^ 
at an apparent magnitude limit Iab < 24. Compa- 
rable figures are achieved by other observ ational cam- 
paigns such as zCO SMOS (ILillv et al.ll200l and DEEP2 
(jPavis et al.l I2003D . being able of surveying selected 
galaxies in various photometric bands ([/, B, R, I) to an 
apparent magnitude limit of about 24. Going deeper in 
redshift, Lyman break galaxy redshift surveys are being 
successful in efficiently getting high quality spectra of 
hundreds of galaxies in t he redshift range 2 .5 < z < 3.5 
within a FOV - Ideg^ (jBielbv et al.ll20Tol) . To get an 
idea, the VIMOS spectrograph can take ^ 500 high qual- 
ity spectra per pointing with an integration time of about 
4h, within a 7x8 arcmin^ FOV, which is coincidentally 
of the same order of the typical error box for a z = 1 
GW event. The typical redshift accuracy of the spec- 
tra is Az < 10-3 (3 X 10-4 in the zCOSMOS survey, 
2 X 10-'^ in the Lyman break galaxy survey), well below 
the typical redshift scale we are interested in (z ^ 0.01). 

Such figures witness the feasibility of efficient spectro- 
scopic redshift determination of a large sample of galaxies 
at faint apparent magnitude {rrir ~ 24), as required by 
our problem. Future spectroscopic survey as BigBOSS 
(jSchlegel et al.l I2009D are expected to further improve 
such figures of merit; a new spectrograph will be able 
to simultaneously get up to 4000 spectra within a single 
pointing of a 7deg FOV. Getting few thousand spectra 
of objects falling within the GW error box in the redshift 
range of interest may be possible in a single observing 
night. At a rrir = 24 cut-off magnitude we generally have 



few hundred to few thousands galaxies in the GW error 
box, but we go deeper (i.e., rrir = 26, feasible with future 
surveys) , the number of spectra may increase drastically. 
For some of the error boxes, we count up to 10"" galaxies 
with nir < 26. However, the requirement of a factor of 
ten more spectra, does not correspond to a significant 
improvement of the results. This is a consequence of 
the self similarity of the galaxy distribution: as long as 
there are enough galaxies in the error box to recover the 
clustering information, the results are basically indepen- 
dent on the assumed cut-off magnitude. A survey with 
a cut-off magnitude of rrir = 24 may indeed be a good 
compromise between reliability of the results and time 
optimization in terms of follow-up spectroscopy. 

The magnitude cut-off defines the number of neigh- 
bor observable galaxies. This is the only practical way 
to weight each galaxy with a local density, rigai = rim^ 
(the subscript rrir refers to the adopted magnitude limit) 
along the lines discussed in Section . Once we have a 
spectroscopic galaxy sample, each galaxy in the error box 
comes with the prior probability to be the host propor- 
tional to rirn^, so the astrophysical prior in equation ([5]) 
could be written as 

p(fi, z) = ^ rim^^i (5(17 - 17j)(5(z - Zi) (11) 

i 

where the sum is over all observable galaxies in the error 
box and fl is the geodesic distance on the celestial sphere 
from the center of the box. At redshifts z > 1 the prior 
probability p(f2, z) becomes almost a continuous function 
(as the example in figure |4]). 

4.3. Approximations and caveats 

Before jumping to the results, we want to mention 
some corrections we made to accommodate the limita- 
tions of our simulations. Firstly, we interpreted one of 
the directions in the snapshot (along the side of the cylin- 
der) as distance from the observer. This is a good ap- 
proximation only if the error box size is small. For large 
error boxes, a uniform distribution in the comoving dis- 
tances does not translate into a uniform distribution in 
redshifts: there is an artificial slope with a bias toward 
low values of z. We have corrected for this slope. Sec- 
ondly, the dumpiness evolves with redshift, which is not 
the case if we use a single snapshot and interpret one 
of the directions as a redshift. To properly account for 
this, we should glue snapshots together and perform an 
interpolation between them. However we wanted to sim- 
plify the setup for this very first attempt. The main idea 
was to check whether the density contrast within the er- 
ror boxes is sufficient to constrain further the error on 
w. If the distribution of density within the error box 
is uniform then we do not gain any useful information. 
However there is a natural bias: for a given measurement 
of Dl, the galaxy further away (larger z) constrains w 
better than galaxy at lower redshift. One can see it from 
the fact that deviation between the curves in D^^z plane 
corresponding to the small deviation in w is bigger for 
large z. This could be counterbalanced by the decreas- 
ing density contrast at large redshift. Here, we corrected 
the slope of the posterior Pj{w\s) by demanding that a 
uniform distribution Pj{9, (j), z) returns a posterior on w 
equal to the prior, i.e., Pj(w\s) = po(w). 
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5. CONSTRAINTS ON THE DARK ENERGY EQUATION OF 

STATE 

In this section we present the results of our simulations. 
Wc tried several setup of the experiment by using differ- 
ent thresholds on the observable apparent magnitude of 
galaxies, different prescriptions for the measurement er- 
rors, and different cosmological priors. For each setup, 
we performed either 100 or 20 realizations of the MBH 
binary population as observed by LISA, together with 
the follow up spectroscopic survey of the galaxies in all 
the error boxes. 



5.1. Fiducial case 

We consider in this subsection 100 realizations which 
we refer to as our fiducial case. For this setup, we 
limit spectroscopic identification of galaxies in the er- 
ror box to an apparent magnitude of to^ £ 24, the er- 
rors in sky localization and in the luminosity distance 
are estimated according to the inspiral part of GW sig- 
nal only, and t he we ak lensing uncertainty is taken from 
IShapiro et all ([20Tnl ). The prior po{w) was assumed to 
be uniform [/[— 0.3 : 0.3] with an exponential decay at 
the boundaries. Such interval is consistent with cur- 
rent 2(7 (95% confidence level) con straints on w {w — 
-0.12 ± 0.27. iKomatsu et al.ll20Tol ). obtained by cross 
correlating seven-year WMAP data with priors com- 
ing from independent m easurements of Ho an d barionic 
acoustic oscillations (sce lKomatsu et aI1l2010l and refer- 
ences therein for full details) , under our same assumption 
for the dark energy equation of state, uj = — 1 — w, where 
w is a constant. Such range is reduced by a factor of 
almos t three (w = —0.0 2 ±0.1) when type la supernovae 
data (jRiess et al.lll998l ) are included. Here we show that 
GW measurements offer a competitive alternative to type 
la supernovae, placing an independent constraint on the 
dark energy equation of state. 

We find that in almost all cases we improve the con- 
straints on w, in other words, the posterior distribution 
is narrower than the prior. Few events at low redshift 
usually play a major role in the final result. One typi- 
cal realization is plotted in the top panel figure \E\ Wc 
split the contribution to the posterior distribution P{w) 
in redshift bands: z £ [0 : 1] (second plot from the left), 
[1 : 2] (third plot), [2 : 3] (fourth plot). Their relative 
contribution and the resulting posterior (black) is given 
in the leftmost plot. In this example the final poste- 
rior probability is almost completely determined by few 
events at low redshift. The second realization, shown 
in the lower panels of figure [51 demonstrates how low 
redshift contributions could give inconclusive results. In 
this particular case, there are two maxima with prefer- 
ence given to the wrong one. The contribution from high 
redshift events could change this ratio as it is shown in 
this example. In many cases the mergers above redshift 
z ~ 1 can constrain w only to a 0.1-0.15 accuracy, but 
they almost always add up coherently giving a maximum 
at the right value (w = 0). This usually helps in case the 
low redshift events return a multimodal P(i«), and is, in 
turn, the power of our statistical method. 

We characterize the results of each setup (100 or 20 
realizations) using the figures of merit shown in figures [HI 
and [71 The first one (figure [6]) is obtained by adding the 
posterior distributions P{w) of all the realizations. We 




0.3-0.15 0.15 0.3 



Fig. 5. — Posterior distribution for w for two particular realiza- 
tions (top and bottom row). In each row, the left plot shows the 
full posterior from all GW events (black curve) as well as contribu- 
tions from different redshift bands. The three right plots show the 
individual contribution for the three redshift ranges, as labelled in 
the panels. 

fit the resulting curve with a Gaussian, characterizing 
the result using its mean wq and standard deviation a^- 
The second figure of merit (figure [7]) shows the result 
of Gaussian fits performed on each individual realization 
(vertical index i): the mean wo{i) is shown as a circle 
and the standard deviation (T^{i) is the error bar. The 
first figure of merit gives collective information, showing 
how well, on average, an individual realization can be 
approximated by a Gaussian fit, while the second figure 
of merit shows the dispersion of the posterior distribution 
across the individual realizations. 

The fiducial case, featuring 100 realizations, is shown 
in panel (a) of both figures [51 and [71 The parameters 
of the global fitting Gaussian mean are wq = 0.0008 and 
~ 0.036, corresponding to a factor of four improve- 
ment in the estimation of w with respect to our standard 
2(7 [—0.3 : 0.3] prior. However the distribution has clearly 
some outliers, recognizable as non-Gaussian tails in fig- 
ure [HI and pinned down in figure [71 For the fiducial case, 
84% of the realizations have a mean value close to the 
true one, i.e. |?«o(*) ~^truc| < 0.1 with an appreciable re- 
duction of the prior range, i.e. (7u,(i) < 0.15 {i = 1, .., 100 
is the realization index). Moreover, most of the outliers 
can be corrected as we will explain in Section [5.61 

5.2. Removing "electromagnetic counterparts" 

Our goal is to demonstrate that we are able to con- 
strain the dark energy equation of state without directly 
observing electromagnetic counterparts. However, for 
some of the low redshift events, the error box is so small 
that only one or two galaxies fall within it. Having one 
or two galaxies in the error box essentially implies an 
electromagnetic identification of the host, so we decided 
to re-analyze the fiducial case removing all such fortu- 
nate events (usually 0-2 in each realization). The fiducial 
case without clearly identifiable hosts is presented in the 
panel (b) of figure [6l Clearly, our results remain almost 
unchanged, the posterior distribution is slightly wider 
(larger sigma) and non-Gaussianity is more pronounced. 

5.3. Choice of the prior for w 
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(a) fiducial case 



(b) without electromagnetic counterpart 



(c) with offset gaussian prior 




(d) apparent magnitude <26 
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Fig. 6. — Collective figures of merit of our experiment. In each panel, corresponding to a diff'erent setup of our experiment as labelled in 
figure, the red solid curve corresponds to the data, i.e. the sum of the posterior distributions of to over all realizations. The blue dashed 
curve is a Gaussian fit with parameter given in the legend of each plot. 
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Fig. 7. — Mean values and standard deviations resulting from 
the Gaussian fit of the posterior P{w). The setup of each panel 
correspond to the one adopted in the panel of figure [6] labelled by 
the same letter). 

Here and in the next subsections we make use of 20 
selected realizations, which we found to be sufficient to 
depict the relevant trends of the analysis. We took f5 
"good" (mean values close to the true and small rms 
errors) and 5 "bad" cases from the fiducial setup. 

In this subsection we study the effect of the prior po{w) 
on the posterior distribution. We considered an extreme 
case: a Gaussian /^{wq = —0.2, a — 0.3). As shown in 
panel (c) of figure [6l the global posterior distribution is 
still centered at the true value w = 0. This demonstrates 
that the final conclusion is basically unaffected by the 
choice of the prior (as long as the prior covers the true 
value) and GW observations, in principle, could be used 
as an independent mean of estimating w. 

5.4. Using deeper surveys 

Here we study the dependence of our results on the 
depth of the follow up spectroscopic survey: i.e. on the 



observability threshold. We considered the same 20 real- 
izations as in the previous section, but now with different 
limits on the apparent magnitude of observable galaxies: 
rrir = 24, 25, 26. The case = 26 is given in panel (d) 
of figures |6] and [T] The results are comparable to the 
fiducial case. They show a small improvement in sigma 
and slightly larger bias for the combined distribution. 
We also notice that 4 out of 5 "bad" cases remain bad. 

We should say few words about the number of galaxies 
used here. As mentioned above, the typical number of 
galaxies for the fiducial case {rur = 24) is less than few 
thousand for events at z < 1 and less than few tens of 
thousands for the high redshift event. For the improved 
observational limit (m^ = 26), these numbers are 2 to 10 
times larger. The fact that our results are not sensitive to 
the depth of the survey reflects the self-similarity of the 
spatial distribution of galaxies in different mass ranges. 

5.5. Improving the sky localization and the luminosity 
distance estimation 

In our fiducial setup, the assumed source sky localiza- 
tion and luminosity distance error are rather conserva- 
tive. In this subsection we consider the effect of improv- 
ing such measurements. So far, we considered only the 
inspiral part of the GW signal; the inclusion of merger 
and ringdown will improve the localization of the source 
by at least a factor of two (^Ic Williams et al. 2010), due 
to the large gain in SNR. We artificially reduced the sky 
localization error coming from the inspiral by a factor of 
two (factor of four in the area), assuming that this will 
be the case if we take the full GW signal. We reanalyzed 
the same 20 realizations with this new error on the sky. 
Because the size of the error box is smaller, the number 
of potential counterparts is reduced by a factor of 4 
compared to the fiducial case. The results are presented 
in panel (e) of figure H) We see that the main effect of 
a better GW source localization is to reduce the number 
of outliers and to remove the non- Gaussian tails in the 
combined probability. As it is clear form panel (e) of 
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figure [71 the main gain comes from improvement of the 
"bad" cases. 

We now consider another estimation of the mean weak 
lensing contribution to th e luminosity distance error, 
given in IWang et alJ (j2002| ) (green square-line curve on 
figure [2]). We take this in combination with improved 
source localization on the sky coming from taking into 
account the merger (as discussed above). We consider 
the same 20 realizations. Results are shown in panel 
(f) of both figures [6] and [T] The improvement with 
respect to all the other cases is obvious. Because the 
marginalized likelihood TTj coming with each galaxy is 
narrower due to the smaller error in the luminosity dis- 
tance, the final posterior on Pj{w) is also narrower. The 
standard deviation (T„, is improved by more than 40% as 
compared to the fiducial case. The non-Gaussian tails 
have almost completely disappeared, due to the removal 
of the outliers (further improvement of the "bad" cases, 
the remaining bad case will be treated in the subsection 
15.61 see also the top panel of figure [8]). With this model 
of the mean weak lensing contribution and assuming the 
full GW signal, the estimation of w is improved by a 
factor of ~ 8 as compared to the initial uniform prior. 

5.6. Consistency check. 

As we mentioned above, some nearby GW event could 
seriously bias the final posterior. We also mentioned that 
the odds for the host to be in a low density region of the 
Universe are not small. The posterior probability P{w) 
reflects the distribution of the mass defined by the astro- 
physical prior Pj{9, 4>,z). A nearby GW event hosted in 
the low density environment could seriously damage the 
final result. An example is given in the top left panel of 
figure [H In order to eliminate or at least test such un- 




FlG. 8. — Each row of panels show our self-similarity eheck for a 
selected realization. In each row, the solid curve on the left panel 
corresponds the final posterior P(w) The solid curves on the right 
panel are the posteriors after removing one event, Pj^{w). 

fortunate cases we performed a self-consistency test on 
our results. Basically we remove one GW event from the 
analysis and see if the resulting posterior P^. (w) distribu- 
tions are consistent. We defined the posterior of all the 
events minus one as: 



Pkiw) 



Jpo{w)Ylj^kPAs\w)dw' 



(12) 



If Pk{w) give similar results for all fc, then we can be 
confident that the result is not biased by one particular 
unfortunate event, and this increases our trust in the fi- 
nal posterior distribution. If, conversely, all Pk{w) but 
one are consistent, then we say that this one event is not 
in line with the remaining events and should be aban- 
doned. In the top panels of figure [8] we see that remov- 
ing one event at low redshift changes the final probability 
completely; the solid (red) line in the right panel is the 
new posterior distribution, consistent with the true value 
w = 0. However, there are still few cases where the self- 
consistency test is not conclusive, and one of them is 
shown in the lower panels of figure HI In this case, re- 
moving one "bad" nearby event produces the red curve 
centered at w = 0, but removing another ("good") event 
results in the green curve, which are mutually not consis- 
tent at all. Since in real life we will not know which event 
is "good" and which one is "bad" , we will not be able to 
make a clear definite statement, and our answer will be 
bi-modal with a probability attached to each mode. 

5.7. Comparison with the optimal case: detection of 
electromagnetic counterparts. 

For comparison, we have also considered the best pos- 
sible case, in which the redshifts of the GW source hosts 
are determined unambiguously through the identification 
of a distinctive electromagnetic counterpart. In this case, 
the redshift of each GW event is known exactly (within 
negligible measurement errors). Therefore, the error on 
w comes only from the error on luminosity distance (GW 
error measurement plus weak lensing). Considering 20 
realizations with a configuration equivalent to the fidu- 
cial case (Section 15. 1[) . the global posterior distribution 
is a Gaussian centered at = with a^j = 0.021 (for 
comparison, see panel (a) of figure [H]). With a configu- 
ration equivalent to our improved case, i.e. better weak 
lensing (Section l5.5p . we obtain = 0.012 (for com- 
parison, see to panel (f ) of figure [6]) . In both case the 
difference between our statistical method and the best 
possible case (all electromagnetic counterparts detected) 
is only about a factor 2. 

6. SUMMARY 

In this paper, we presented a statistical method for 
constraining cosmological parameters using LISA obser- 
vations of spinning massive black hole binaries and red- 
shift surveys of galaxies. Our approach does not require 
any direct electromagnetic counterpart; instead, the con- 
sistency between few dozen of GW events imposes con- 
straints on the redshift-luminosity distance relationship. 
This, in turn, allows us to estimate cosmological param- 
eters. This method strongly relies on the non-uniformity 
(i.e., clustering) of the galaxy distribution within the un- 
certainty error box set by LISA observations, weak lens- 
ing and priors on the cosmological parameters. 

For this first exploratory study, we fixed all the cosmo- 
logical parameters but one, it;, describing the effective 
equation of state for the dark energy. We used the Mil- 
lennium simulation to model the Universe at different 
redshifts. We used a particular (VB) hierarchical MBH 
formation model to mimic the MBH binary population 
observed by LISA. Using this setup, we considered be- 
tween 20 and 100 realizations of the observed LISA bi- 
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nary population. We tried two different models for es- 
timating the error in luminosity distance due to weak 
lensing, we also looked at the effect of including merger 
and ringdown via improvement of the sky localization. 
Wc checked the robustness of our final result against dif- 
ferent depth of future spectroscopic galaxy surveys. 

Our fiducial case, based on conservative assumptions, 
shows that we are able to constrain w to a 8% level {2a), 
i.e., we improve its estimate by a factor of 4 as com- 
pared to the current 95% confidence interval obtained 
by cross correlating the seven-year WMAP data analysis 
with priors coming fr om Hp measurements and barionic 
acoustic oscillations (jKomatsu et al.l [20101 ). Such new 
measurement would be at the same level (25% better on 
average) than current constraints based on seven-year 
WMAP data plus type la supernovae observations. The 
optimistic case (smaller weak lensing disturbance and full 
GW waveform) allows us a further improvement by an- 
other factor of two, providing a factor of ~ 2.5 tighter 
constraint than current estimates including supernovae 
data. Our results are most sensitive to the weak lensing 
error (witnessing once more how critical is the issue of 
weak lensing mitigation for cosmological parameter esti- 
mation through GW observations) and are almost inde- 
pendent on the depth of the rcdshift survey (provided we 
have a reasonable number of rcdshift measurements per 
error box). 

In the majority of the realizations the most information 
comes from few events at low redshift, and high rcdshift 
events do help in case of multimodal structures in the 
posterior distribution. We suggested a self-consistency 



check based on the similarity of the posterior distribution 
from each GW event. This increases our confidence in 
the final result and allows to reduce the risk of incurring 
in unfortunate outlier realizations for which we can not 
place useful constraints on w. We also compared our sta- 
tistical method to the optimal situation in which electro- 
magnetic counterparts to the GW sources are identified, 
finding an improvement of a factor of two in the latter 
case. In absence of distinctive electromagnetic counter- 
parts, statistical methods like the one presented here can 
still efficiently constrain cosmological parameters. 

Although the main result of the present paper is en- 
couraging, it was obtained assuming a fixed cosmological 
model with one free parameter only: the w parameter de- 
scribing the dark energy equation of state. Even though 
we will likely have a good knowledge of most of the other 
cosmological parameters by the time LISA will fly, it is 
worth considering models with more degrees of freedom. 
In following studies, we intend to consider a more realis- 
tic situation by releasing other cosmological parameters, 
testing LISA capabilities of setting constraints on a multi 
parameter model. 
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